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ABSTRACT 

Dust at the midplane of a circumstellar disk can become gravitationally unstable and fragment 
into planetesimals if the local dust-to-gas ratio fio = Pd/Pg is sufficiently high. We simulate how 
dust settles in passive disks and ask how high //o can become. We implement a hybrid scheme that 
alternates between a ID code to settle dust and a 3D shearing box code to test for dynamical stability. 
This scheme allows us to explore the behavior of small particles having short but non-zero stopping 
times in gas: < i s top *C the orbital period. The streaming instability is thereby filtered out. Dust 
settles until Kelvin-Helmholtz-type instabilities at the top and bottom faces of the dust layer threaten 
to overturn the entire layer. In this state of marginal stability, fio = 2.9 for a disk whose bulk (height- 
integrated) metallicity Ed/Eg is solar — thus increases by more than two orders of magnitude from 
its well-mixed initial value of /L«o,mit = Ed/£ g = 0.015. For a disk whose bulk metallicity is 4x solar 
(Ato,init = 5]d/S g = 0.06), the marginally stable state has /iq = 26.4. These maximum values of fio, 
which depend on the background radial pressure gradient, are so large that gravitational instability of 
small particles is viable in disks whose bulk metallicities are just a few (<; 4) times solar. Our result 
supports earlier studies that assumed that dust settles until the Richardson number Ri is spatially 
constant. Our simulations are free of this assumption but provide evidence for it within the boundaries 
of the dust layer, with the proviso that Ri increases with Ed/£ g in the same way that we found in 
Paper I. Because increasing the dust content decreases the vertical shear and increases stability, the 
midplane fj,Q increases with Ed/£ g in a faster than linear way, so fast that modest enhancements in 
Ed/E g can spawn planetesimals directly from small particles. 

Subject headings: hydrodynamics — instabilities — planets and satellites: formation — protoplanetary 
disks 
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1. INTRODUCTION 

Dust can settle quickly in gaseous protoplanetary 
disks. In a passive (non-turbulent) nebula, a particle's 
vertical height z above the midplane obeys 



(i) 



where the first term on the right-hand side accounts for 
gas drag, and the second term accounts for stellar grav- 
ity when z <C r, the cylindrical radius. Here fix is the 
Keplerian orbital frequency and 



stop 



W« rc i 



(2) 



is the momentum stopping time of a particle of mass m 
moving at speed v le \ relative to gas. Express i ons fo r the 
drag force Fp can be found in lAdachi et al.l (|1976l ) and 
iWeidensch illing (19771). We are interested in small, well- 
coupled particles having non-zero stopping times much 
shorter than the dynamical time: < r s = f^K^stop *C L 
Spherical particles of radius s and internal density p s that 
experience Epstein drag (Frj cx s 2 v IC \ so that t stop does 
not depend on i> re i) settle to the midplane at terminal 
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For this and all other numerical evaluations in this pa- 
per, we use a background disk that is F times more 
massive than the minimum-mass nebula of Chiang & 
Youdin (2010, hereafter CY10; see Appendix [SJ. For 
such a disk Nettie is nearly independent of stellocentric 
distance. The assumption that particles are spherical 
may not be too bad because fractal aggregates of grains 
are ex pected to compactify as they collide with one an- 
other (IDominik fc : Tielensl Il997t IDuUemond fc Dominik! 
[20051 [Qrmel et al.ll2007l) 

For millimeter-sized particles, the settling time £ se ttie 
is much shorter than the disk lifetime, me asured in Myr 
(|Hillenbrandl [2001 iHernandez et all [20081) . By compar- 
ison, micron-sized and smaller particles stay suspended 
at least one scale height above the midplane as long as 
the gas disk is present. To the extent that collections 
of particles of different sizes tend to place their mass 
at the upper end of the size distribution and their sur- 
face area at the lower end, we can expect most of the 
solid mass in disks to sediment out into a thin sublayer, 
leaving behind the smallest of grains to absorb incident 
starlight in a flared disk atmosphere. On the whole this 
picture is consistent with observed spectral energy dis- 
tributions of T Tauri disks, although some models hint 
that large grains might remain lofted up in a disk two gas 
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scale heights thick (D'Alcs sio et al.l [20061. Settling can 
only proceed when and where disk turbulence dies, in 
regions where gas is insufficie ntly dense to sustain grav- 
itoturbulence (|Gammiell2001h and too poorly ion ized to 
be magnetorotationally unstable (jGammiel ll 996) . Our 
current understandin g of disk turbulence easily admits 
such passive regions. iTurner et all pOlOh found m nu- 
merical simulations that even when disk surface layers 
were magnetorotationally unstable, grains at the mid- 
plane set tled much as they would in a laminar flow. 
Recently iPerez-Becker fc Chiang! (|2010l ) estimated that 
practically the entire disk would be immune to the mag- 
netorotational instability because ion densities are too 
low for the plasma to drive turbulence in the overwhelm- 
ingly neutral gas. 

In passive disk regions, how far down do dust par- 
ticles settle? When particles are small enough not to 
be affected by aero dynamic streaming instabilities (e.g., 
iBai fc Stone! l2010"ah . we expect them to settle until the 
dust density gradient dpa/dz becomes so large, and 
the consequent vertical shear in orbital velocity dv^/dz 
so strong, that the sublayer is on the verge of over- 
turning by a Kelvin -Hclmholtz-type instability (KHI; 
Wcidcnschilling 1980). An order-of- magnitude estimate 
of the minimum layer thickness can be derived using the 
Richardson number 



Ri 
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which if less than some critical v alue Ri r . r \t may signa l 
that the layer is KH unstable (e.g,. IDrazin fc R cid 2004). 
Here g is the vertical gravitational acceleration and p = 
Pd + Pg is the total density of dust plus gas. In Lee 
et al. (2010, hereafter Paper I), we found that Ri C rit 
increases with Ed/E g , the ratio of dust to gas surface 
densities, a.k.a. the bulk metallicity. For disks of bulk 
solar metallicity, we determined empirically that _Ri crit ss 
0.2. 

To translate the Richardson number (j4j into a critical 
dust layer thickness, first recognize that the orbital veloc- 
ity Vcj, depends on the local dust-to-gas ratio p = Pd/Pg 
according to 



= tt K r 1 



in the inertial frame, where 
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is a dimensionless measure of the strength of the back- 
ground radial pressure gradient dP/dr , with gas scale 
height H g and sound speed c s (e.g., INakagawa et al.l 
1986). When dP/dr < 0, pressure provides extra sup- 
port against radial stellar gravity and so drives the gas 
to move on slower than Keplerian orbits. The orbital 
velocity depends on p as in ([5]) because dust-laden gas, 
weighed down by the extra inertia of solids, is accelerated 



less by the radial pressure gradient than is dust-free gas, 
and so must hew more closely to Keplerian rotation. Call 
the critical layer height Azr,; for which Ri = Ri CI n, and 
assume the midplane gas-to-dust ratio po ^> I (above 
the layer p <C 1). Then evaluating equation Q with 
the approximations g ~ — Sl^AzRi (no self- gravity), 
p~ 1 dp/dz ~ — l/Azpu, and dv^/dz ~ — y^Kx/AzRi, we 
find 
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Equation (0 indicates the dust layer could be quite 
thin, subtending on the order of 1% of the gas scale 
height. Is this thin enough for the dust to self-gravitate 
and hopefully fragment into planetesimals? One can 
compare the midplane density to the "Toomre den- 
sity" required for the disk to undergo gravitational in- 
stability on the dynamical time fl^ 1 (jSafronovl 119691 : 
IGoldreich fc Wardl [1971 Fl As reviewed by CY10, the 
Toomre density it0 
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where the numerical evaluation is for a central stellar 
mass M* equal to 1A/ Q . Now the actual midplane (sub- 
script 0) density is 



Po = Pgo+Pdo = 2.7x 10- 9 F(l+ M o) (-^j ) 



r \ -39/14 



g cm 



which means the midplane dust-to-gas ratio must be 

-3/14 
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for the midplane density to match the Toomre density. 
By comparison, in our crude model of a dust sublayer 
whose height above the midplane cannot be smaller than 
AzRi, the midplane dust-to-gas ratio cannot exceed 



E d /(2Az Ri ) 
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E d /E g \ / 0.2 
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p g o V 0.015 J \Rimt, 

which is nominally smaller than p ^ Toomre by more than 
an order of magnitude. Here the bulk (height-integrated) 
metallicity Ed/E g is normalized to solar abundance 

5 Even if the Toomre density is attained so that marginally 
bound clumps of dust and gas form in time, continued col- 
lapse of dust is not guaranteed. For dust to concentrate further 
it must sediment to the centers of the gas clumps over timescales 
^settle 3> • D urm g this time the dust clumps are susceptible to 
erosion by gas streaming or turbulence (e.g., Cuzzi ct al. 2008). 

6 Strictly speaking, the Toomre criterion for gravitational 
instability is derived for two-dimensional disk s charac t erized 
by surface densities, not v olume densities ijToomrel 119641 ; 
Goldrcich & Lyndcn-Bcll 1965). To derive our Toomre volume 
density, we assign a half-thickness to the disk equal to c/Hk, where 
c is the velocity dispersion of the dust + gas mixture. This assign- 
ment is not rigorous; see CY10. 
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(|Loddersll2003l) . assuming all metals have condensed into 
grains 

For many years the fact that /io.Hi falls short of 
A*o,Toomrc was believed to rule out the formation of plan- 
etesimals by collective effects, self-gravitat ional or oth- 
erwise (e.g.. iWeiden schillin g fc Cuzzilll993T) . But there 
are more ways to achieve the Toomre density than verti- 
cal settling. A dissipative form of gravitational instabil- 
ity can, in principle, collect particles radially into over- 
dense rings eve n when self-gravity is weaker than stel- 
lar ti dal forces (lWardl ll976t IWardll2000t iCoradini etlll 
119811 ; lYoud in 2005; for a sim ple explanation of the insta - 
bility, see the introduction of lGoodman fe Pindori 2000). 
It is not clear whether this instability, which operates 
over lengthscales and timescales longer than those char- 
acterizing the Toomre instability by at least a factor of 
(MOjToomro/Mo) 2 ! can compete wi th other effec ts that seek 
to rearrange dust and gas (e.g., Youdin 2005). 

Another alternative is to invoke larger dust parti- 
cles that are only marginally coupled to gas; these 
can clump by the aerodynam ic streaming instability 
(SI: lYoudin fe Goodman! 120051: Uohansen et all [20091: 
IBai fe Stond l2010at IBai fe Stondl2010cri. I n their3D nu- 
merical simulations. Uohansen et al.l (2009) reported that 
particles having t s = 0.1-0.4 — corresponding to sizes of 
a few centimeters at r = 5 AU if F = 1, and larger 
sizes if F > 1 — concentrated so strongly by aerodynamic 
effects that planetesimals effectively hundreds of kilome- 
ters across coalesced within just a few orbits. To obtain 
this result, Uohansen et al.l (|2009l ) initialized their simu- 
lations by placing the bulk of the disk's so lid mass into 
particle s approaching decimeters in size. IBai fe Stone! 
(2010a) greatly expanded the range of r s modeled and 
found similar results for their 3D simulations: in the 
highly turbulent states driven by the SI, instantaneous 
densities exceeded the Roche densitj0 when the disk's 
solids were all composed of particles having r s = 0.1- 
1 and the bulk height-integrated metallicity was about 
twice solar; see run R10Z3-3D in their Figure 5. For 
this same run, the time-averaged dust-to-gas ratio at the 
midplane was ~12, a factor of a few less than the Toomre 
threshold; see their Figure 4 and compare with our equa- 
tion (flU)) . By contrast, when half or more of the disk's 
solid mass had r s < 0.1, or when disks had smaller metal- 
licities, their simulated densities fell short of the Roche 
and Toomre densities by more than an order of magni- 
tu de. Note that we are quoting from the 3D simulations 
of IBai fc Stoiii (I2010afl . 

Given how sensitive the SI is to the existence of 
marginally coupled particles (centimeter to meter sized 
for t s ~ 0.1-1, r ~ 1-30 AU, and order unity F), 
whether enough such particles actually exist in proto- 
planetary disks for the SI to play a dominant role in 

7 The assumption that all metals have condensed is valid only 
for temperatures T < 41 K. For 41 K < T < 182 K, methane 
ice sublimates and £ d / s g « 0.78 x 0.015 (Loddcrs 2003; CY10). 
For our fiducial disk described in Appendix IAI T < 182 K for 
r > 0.4 AU. For T > 182 K, water and other ices sublimate and 
the maximum E^/Sg RJ 0.33 X 0.015. Reductions in E,j/Eg due to 
sublimation may be offset by radial pileups of dust. 

8 The Roche density is that required for a fluid satellite to be 
gravitationally bound against tidal forces exerted by a central body. 
It is greater than the Toomre density by a factor of ~7ir, and is the 
more appropriate threshold density for the highly localized clumps 
of dust created by the SI. 



planetesimal formation remains an open and delicate is- 
sue. Appeal is often made to observed spectral energy 
distributions and images of T Tauri disks at centime- 
ter wavelengths; these suggest that much of the solid 
mass is in millimeter to centimeter sized particles (e.g., 
ID'Alessio et"al|[200lt IWilner et al.l [2005h . Larger sized 
particles are plausibly also present but are not inferred 
for want of data probing the disk at longer wavelengths. 
One problem concerns how quickly t s <; 0.1 particles can 
be grown, and how they can su rvive orbital decay b y gas 
drag. In the 3D simulations of IBai fc Stond ()2010aD . the 
SI clumped particles strongly enough for self-gravity to 
be significant when t s ;> 0.1 particles comprised more 
than half of the disk's solid mass. It is unclear whether 
particle-particle sticking can build up such a population 
before it is lost to the star by gas drag. This concern 
is ameliorated by enhancements in particle density (pile- 
ups) that may occur as particles drift radial ly inward 
(|Youdin fc Shul l200l lYoudin fc Cmangl l200l . and by 
the reduction of drift speeds brou ght about by multiple 
particle sizes (|Bai fc Stondl2010at see their Figure 8). 

Regardless of which scenario nature prefers — particle 
concentration by the streaming instability; dissipative 
gravitation into rings; or dynamical collapse of a verti- 
cally settled sublayer, which is the subject of this paper — 
all the proposed ways of forming planetesimals depend 
on knowing how far down dust settles and what maxi- 
mum dust-to-gas ratios /io can be attained at the mid- 
plane. Our order-of-magnitude estimate in equation (jlll) 
requires testing. Among the mo st realistic simu l ations 
of p article settlin g are th ose by Uohansen et al.l (|2009l ) 
and IBai fc S tond (I2010af). both of which concentrated 
on the SI. Uohansen et alj (|2009f) reported that super- 
centimeter sized particles settled into sublayers in which 
the midplane-averaged /io ranged from 0.6 to 9.0 as the 
bulk m etallicity ranged from 1 to 3x solar. IBai fc Stond 
(|2010af) found that the highest r s particles settled the 
most, driving turbulence that lofted smaller r s particles 
to greater heights. They argued that in their simulations, 
all of which were characterized by maxr s > 0.1, particles 
were so strongly stirred by the SI that the KHI never 
manifested. 

To complement these studies, we would like to under- 
stand the settled equilibrium states of disks composed 
entirely of small particles, well but not perfectly coupled 
to gas (0 < maxT s < 1), isolated from the complicating 
effects of the streaming instability but not other instabil- 
ities like the KHI. Some previous attempts in this regard 
relied on assumed forms fo r the density profile of set- 
tled dust. Barranco (2009) presumed the dust density 
profile was Gaussian in shape, and did not seek to de - 
termine the maximum value o f jip pe r se. Chiand (120081) 
and P aper I, following [Siki^ <1 1998ft and lYoudin fc Shul 



( 2002), assumed the dust density profile had a spatially 
constant Richardson number. We found in Paper I that 
under this assumption, in a disk of bulk solar metallicity 
(E(j/Sg = 0.015), the sublayer could remain KH stable 
for /io as high as 8 — a value that is nearly an order of 
magnitude higher than our crude estimate in (|11[) . and 
as such lowers the hurdle to forming planetesimals by 
g ravitational instability. 

IWeidenschilhnd" (pOM ) and IWeidenschilling! ((2010ft did 
not assume a shape for the marginally stable density pro- 
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file to which dust relaxes. Instead the profile was calcu- 
lated from a one-dimensional model that balanced the 
downward flux of particles by gravity (including verti- 
cal self-gravity) with the upward flux due to turbulent 
diffusion. The model used prescriptions for KH turbu- 
lence, with a num ber of parameters c hosen to match the 
3D calculations of iCuzzi et al.l (|1993t ) — which also relied 
on a p rescribed form of the turbulence. iWeidenschillind 
( 2006) found that the density of millimeter-sized parti- 
cles could exceed the Toomre density at r = 3 AU in 
a disk approximately 1.3x as massive as our minimum- 
mass solar nebula, and whose bulk metallicity Ed/^g = 
0.054 3.6 x solar. See his Figure 10, but note that 
3.6 x solar metallicity « 16 x his "nominal" abundance 
of solids, and that his "critical" density is 3x larger than 
our Too mre density. 

Like iWeidenschilhngl (puM [2010h . we seek the 
marginally stable state to which dust settles, in the limit 
that particles are well but not perfectly coupled to gas 
(0 < r s -c 1). We improve upon these earlier studies 
by not prescribing or parameterizing the turbulence, but 
by letting turbulence arise and evolve naturally from our 
3D integrations of the standard fluid equations. In par- 
ticular, our work is free of the popular but untested as- 
sumption that dust settles until the Richardson number 
equals a constant everywhere. We allow dust grains to 
fall until they are stopped by whatever instabilities they 
self-generate. In the calculations presented here we as- 
sume that dust begins well mixed with gas in a Gaussian 
density profile, and then follow the dust into whatever 
non-Gaussian distribution it seeks to relax. Although we 
try only a Gaussian initial profile, our method accommo- 
dates arbitrary initial conditions. 

At the heart of our approach lie two codes. The first 
code is in one dimension (z) and computes the vertical 
drift of dust grains at their terminal velocities. Though 
incapable of deciding whether the density profiles it gen- 
erates are prone to the KHI (or any other instability) , the 
ID code can evolve dust profiles for the entire settling 
time ^settle) which can be arbitrarily long for arbitrarily 
small grains. The task of assessing stability is reserved 
for the second code: the spectral, anelastic, shearing box 
code of lBarrancol (|2009| ) which treats gas and dust in the 
perfectly coupled r s = limit. Though incapable of al- 
lowing dust to settle out of gas, the 3D code accounts for 
the complicated interplay of vertical shearing and rota- 
tional effects to decide whether a given dust layer over- 
turns from the KHI (or some other instability). It tests 
dynamical stability by running for dozens of dynamical 
times tdyn = f^K < Our procedure involves alternating 
between these two codes: allowing dust to settle over 
some fraction of the settling timescale t se ttie using the 
ID code; passing the results of the ID code to the 3D 
code and allowing the dust profile to relax dynamically 
over timescales td yn ] passing the results of the 3D code 
back to the ID code for further sedimentation on the set- 
tling timescale; and so on, back and forth, until the mid- 
plane dust-to-gas ratio stops increasing, at which point 
the marginally stable state is identified. 

In <j2]we describe our method in full. Results are pre- 
sented in extended in 21 and summarized and dis- 
cussed in fJS] 



2. METHOD 

As sketched in SJTJ to find the marginally stable state 
to which small dust grains relax, we alternate between 
two codes: a ID code that tracks how dust drifts toward 
the midplane on the settling ti mescale ts e ttie. a id a 3D 
shearing box code developed bv lBarrancol (|2009f) that al- 
lows dusty gas to stabilize on the dynamical timescale 
tdyn "C t B ettle- The 3D code integrates the anelastic fluid 
equations for perfectly coupled dust and gas using a spec- 
tral method. It includes a background radial pressure 
gradient to drive a vertical shear. De t ails about the 3D 
code are in lBarranco fc Marcus! (|2006D . lBarrancol (|2009D . 
and Paper I. 

Dust and gas are initially well mixed with a spa- 
tially constant density ratio: [pd(z) / p s (z)]mit = Minit = 
constant. We set /x, n it equa l to either solar metallicity 
(Minit = 0.015: lLoddersll2003l : see also footnote [7]) or four 
times solar metallicity (ptinit = 0.06). To determine the 
initial form of the dust density profile Pd{z), we solve the 
equation for vertical hydrostatic equilibrium where gas 
is assumed to be initially isothermal: 



p g + pa dz 



whence 



Pd = MinitPgO exp 



(1 + Minit)^ 2 



2H? 



(12) 



(13) 



for constants p lni t, a characteristic initial height H g = 
c s /f2K, and the midplane gas density p gu . 

Taking pi n a to be constant over all gas scale heights is a 
simplifying but probably unrealistic assumption. Even if 
particles of fixed size s carry the bulk of the disk's solid 
mass, such particles would not likely begin well mixed 
with gas everywhere. Densities at altitude may be too 
low to permit particles of size s to coagulate, and to keep 
such particles aloft once formed. Nevertheless the error 
accrued by our assumption of constant p- m it is small inso- 
far as increasingly small amounts of mass are contained 
at larger altitudes. Moreover, we will find evidence that 
the final marginally stable state to which dust relaxes is 
insensitive to initial conditions f ij5.ip . In any case, we 
will point out in fJ3] which features of the evolving dust 
profile are artifacts of our assumed initial condition. 

Equation (fl~3)) defines the initial dust profile used by 
the ID code, whose grid extends from z = to z — 
3-ffg. The ID Lagrangian code uses particles to track the 
motion of dust mass. Each particle represents the same 
amount of dust mass. Any dust density profile Pd(z) can 
be converted into particle positions and back again. The 
closer particles are spaced, the greater is pd- 

Starting with equation fT3"]) , we proceed as follows: 

1. ID code: Initialize positions of dust particles 
and establish hydrostatic equilibrium for the 
gas. 

Given Pd(z), calculate the positions of ^60,000 par- 
ticles in the ID code. Also determine the hydro- 
static gas density p g (z) at each particle's position 
by solving equation (TT2"j) FI 



Actually the calculation of p g (z) can be neglected to good 
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2. ID code: Settle dust particles by one 
timestep At. 

By equating the vertical gravitational force oc z 
to the Epstein drag force i*b oc p g v rc \ (e.g., 
iWeidenschilling 1977), assign terminal velocities 

v Ic i ex — (14) 
Ps 

to each particle. Advect each particle vertically 
downward by a distance i> ro iAi, where At is chosen 
small enough that particles do not overtake one an- 
other. Note that the coefficient of proportionality 
on the right hand side of (fl4"|) . which depends on 
quantities such as the disk mass parameter F and 
grain properties s and p s , does not affect the shapes 
of the density profiles generated so long as it is the 
same for all particles. 

Bin particle positions and recalculate pd(z). 

3. ID code: Repeat (1) + (2) until the mid- 
plane dust density pd(0) rises by 30%. 

4. ID — » 3D code: Insert results of the ID code 
into the 3D code and run the 3D code for ten 
orbits. 

Let z max be the position of the highest particle 
in the ID code, and set the dimensions of the 
shearing box in the 3D code to be (L r , L^, L z ) = 
(1.455, 2.91, 4)z max , resolved by {N r ,N^,N z ) = 
(32,64,128) gridpointsFJ Initialize the 3D code 
by assigning p(z), as calculated by the ID code, to 
each horizontal gridpoint (r, </>)Q The 3D code ini- 
tializes the remaining variables — velocity, gas den- 
sity, temperature, and enthalpy — to ensure dynam- 
ical equilibrium; see section 2.2 of Paper I. 

The background radial pressure gradient is param- 
eterized by the variable w max : 

1 OP 

7T- = 2fi K Vmax ■ (15) 

Pg or 

We fix u max = 0.025c s for all simulations. Phys- 
ically, f max ~ c£/(Qkt) represents the difference 
in azimuthal velocity between pressure-s upp orted 
dust-free gas and a strictly Keplerian flow Pi 

approximation, as the hydrostatic gas density deviates only slightly 
from a Gaussian throughout the evolution. Even when fi > 1 
near the midplane, Ap g /p g ~ (z/ H g ) 2 fi, which for our parameters 
remains much less than unity. In fact, the gas density is practically 
constant once the dust falls to z < 0.1-ffg. 

10 These choices imply that every z max length in the r and <fr 
directions is resolved by 22 grid points. In the z direction we 
need N z = 128 points to achieve comparable resolution because 
the vertical grid differs fr om t he horizontal grid (see section 2.4 of 
Paper I). As discussed in £14.11 we test the robustness of our results 
to box size by using bigger boxes as the marginally stable state is 
approached. 

1 The transfer of fi(z) from the ID code to the 3D code involved 
some smoothing because the vertical grid for the 3D code is ~10x 
coarser than that of the ID code. We captured all features of the 
ID dust profile Pd( z ) to within ~10% for z < 0.8z max . Fractional 
errors generally increased away from the midplane and were largest 
at 2 m ax where the dust content goes to zero. 

12 This value of « max coincides with the standard value from 
Paper I, although technically the minimum-mass disk model in 



Before running the 3D code, perturb p(r, 4>, z) by 
an amount 

A/i(r, </>, z) = A(r, <f)p{z) [cos(7rz/2z max ) 

+ sin(7rz/2z max )], (16) 

where A(r, (f>) is a random variable constructed in 
Fourier space (see the discussion following equation 
31 of Paper I). Fix the root-mean-square of the 
perturbations to be A lms = {A 2 ) 1 / 2 = 10~ 3 . 
Run the 3D code for ten orbits. 

5. 3D code: Assess stability. Extend simula- 
tions beyond ten orbits as necessary to make 
this assessment. 

Label the dust profile "KH-unstable" if the hori- 
zontally averaged dust-to-gas ratio at the midplane 

(p(z = 0)) as a function of t 

decreases by more than 15%. Otherwise, monitor 
the horizontally averaged vertical kinetic energy at 
the midplane: 

(pv 2 (z = 0)) /2 as a function of t. 

If (pv 2 }/2 monotonically decreases or levels off, la- 
bel the dust profile "KH-stable." If (pv 2 )/2 is 
increasing towards the end of the simulation, ex- 
tend the integration an additional ten orbits and 
re-assess stability. Repeat step (5) as necessary. 

6. If "KH-unstable," stop. 

Identify the last KH-stable dust profile, gener- 
ated in the iteration just previous to that of the 
KH-unstable simulation, as the "marginally stable 
state.'H 

7. If "KH-stable," pass results of the 3D code 
back to the ID code and return to step (1). 

Fit a polynomial p, po i y (z) to the final, horizontally 
averaged dust-to-gas ratio (p{z)) as calculated by 
the 3D code. Adjust the order of the polynomial 
to capture all features of the profile. If one poly- 
nomial is insufficient, use two to create a piece- 
wise function. Convert p po i y (z) to Pd{z) by assum- 
ing the gas profile to be Gaussian (see footnote |HJ) : 
p d (z) = Mpoiy(z) • p g0 exp[-z 2 /(2i7|)]. Using this 
Pd{z), return to step (1) for the next iteration. 

Our method and results apply to any location in a 
disk of any mass (arbitrary r and F), provided our input 
assumptions that self-gravity is negligible and w max /c s = 
0.025 are satisfied. They also apply to any particle size 
to the extent that the disk's solid mass is concentrated 
in particles of a single size (so that z and p g uniquely 
determine w ro i; equation [14]), and to the extent that such 

Appendix [K\ gives D max /c s = 0.036. Both values of n max corre- 
spond to cool disks passiv ely heated by their central stars (e.g., 
Chiang & Goldrcich 1997). Smaller fmax results in thinner and 
denser dust laye rs and thus promotes the formation of planctcs- 
imals. See £|5,2I and Appendix [B] for how our results depend on 

^m ax ■ 

13 This marginally stable state will be superseded by the 
marginally stable state identified under an improved scheme in i|4.2l 
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part icles are undisturbed by streaming instabilities (§TJ 
£15. 3[) . Another way of saying all this is to note that our 
calculations are carried out in dimensionless units. 

3. RESULTS 

To orient the reader, in Figure Q] we show results ob- 
tained from the ID code only. The dust is initially well 
mixed with gas at solar metallicity (jUfrut =0.015). As 
dust settles and the midplane dust-to-gas ratio fj,o in- 
creases, sharp cusps appear at the edges of the dust layer 
where particles pile up vertically. Pileups occur because 
particle fluxes Pd|^roi| oc fi\z\ increase with increasing 
height \z\. This follows from our assumption that fii n n is 
constant, which as noted at the beginning of £]2]may not 
be realistic. 

Unlike us, iGaraud fc Linl (|2004[) did not find vertical 
pileups at the edges of their layer because they chose 
their initial dust profile to have a scale height equal to 
0.1-ffg. Their initial /i profile decreased with \z\ more 
quickly than 1 / 1 z | , and thus did not satisfy the condition 
for pileups. We verified this by inserting their initial 
profile into our ID code. 

The shapes of the settled dust profiles n(z) and their 
relative spacing in time are independent of the dust in- 
ternal density p s , dust particle size s, and the scaling 
parameter F for disk mass. Changing these parameters 
only alters the absolute physical time elapsed (equation 
[3]). Relative time is tracked by the dimensionless param- 
eter / = t/t S ettio, labeled on this and many subsequent 
figures. 

Below we compare these ID-only results to those that 
include the full 3D dynamics. The solar metallicity case 
is described i n §3.11 The metal-rich case (/x; n i t = 0.06) is 
presented in £13.21 

3.1. Solar Metallicity 

Figure [2] traces the evolution of dust that starts well 
mixed with gas at solar metallicity. Plotted are several 
KH-stable curves from the 3D code resulting from step 
(5) of our procedure. For ease of comparison with the 
purely ID results, the relative timestamps in Figure [21 
measured by /, coincide with those in Figure [1] The 
leftmost curve at / = 1.0 represents the marginally sta- 
ble state identified using our standard procedure. This 
state achieves a midplane dust-to-gas ratio of = 2.45, 
about an order of magnitude below the value required 
for gravitational instability fequation [T0|). In §J] we ex- 
tend our procedure to see if we might achieve still higher 
dust-to-gas ratios. 

Comparing Figures [1] and [21 we see that the (possibly 
unrealistic) pileups at the edges of the dust layer do not 
survive in the dynamical 3D code. By / 0.44, the 
pileups are nearly gone. At this point, the vertical extent 
of the dust layer z max has shrunk to ^0.1_ff g , and the only 
pileup present is the one at the midplane. 

The instability that eliminates the pileups at the edges 
of the layer is likely related to the Ray leigh- Taylor in- 
stability (RTI), triggered by heavy fluid lying on top of 
lighter fluid, and we will refer to it henceforth as such. 
The RTI originates locally at the edges of the dust layer. 
By contrast, the midplane is relatively stable (at least 
until the marginally stable state is reached). Another 
way of seeing this is to note that midplane dust-to-gas 
ratios fio in Figures Q] and [2] agree to within 25%. Closer 



examination reveals that those in Figure [21 are consis- 
tently higher. This suggests that the RTI transfers some 
of the dust in the pileups to the midplane. 

Figure [3] confirms this transfer mechanism. The top 
middle panel shows that over the course of a 20-orbit- 
long 3D simulation (iteration #6, occurring at a time 
/ = 0.31, out of a total of 19 iterations), dust is re- 
distributed from the layer's edges to the midplane, rais- 
ing by about 20%. ' Note that the effect of the RTI 
has been to transport dust toward the midplane, not to 
higher altitudes. The RTI is confined to where dust is 
unstably stratified (increasing total density in the direc- 
tion opposite to gravity). 

Compare this behavior with that in the top row of 
Figure HI which documents a later iteration, #16. The 
top middle panel shows that an instability has occurred 
near the edges of the dust layer. Dust is redistributed 
to higher, not lower, altitudes. The midplane is not af- 
fected. The instability at this relatively late stage of set- 
tling is probably driven by the vertical shear associated 
with strong density gradients at the edges of the layer, 
and we will refer to it henceforth as the Kelvin- Helmholtz 
instability (KHI). As a result of the KHI, gradients in 
density and velocity are reduced. 

The marginally stable state identified using our stan- 
dard procedure is displayed in Figure [5] The bottom 
panels show that during the last iteration #19, the usual 
30% increase in the midplane fio (left bottom) results in 
a KH- unstable profile (middle bottom). In the top pan- 
els, we redo iteration #19, this time incrementing /iq by 
only 10% (left top). The resultant profile is KH stable 
(mid dle and right top panels), and has (^o) = 2.45. In 
£14.21 we modify our standard procedure and extend it 
to later times to achieve still higher dust-to-gas ratios in 
stable flows. 

The /i-profiles in FiguresEHSbetray oscillations just in- 
side the edges of the dust layer. We believe these ripples 
are artificial because when each first appears, it spans 
only a few grid points of the 3D code: see the / = 0.054 
profile of Figure [21 which shows two nascent ripples. The 
features probably arise because the truncated Chebyshev 
series used to model the flow in z has too few terms 
to adequatel y capture the steep vertical density gradient 
(Gibbs 1898). Originating in the 3D code, the ripples are 
then amplified as mini-pileups in the ID code. We could 
have tried to smooth away these oscillations by reducing 
the order of our polynomial fit (step 7 of our procedure) , 
but chose instead to retain all features of the dust pro- 
file generated by both codes to minimize bias. In any 
case the oscillations are eventually erased by instabilities 
during the later stages of settling (Figure [2]). In and of 
themselves the oscillations do not appear to introduce in- 
stabilities, which as discussed above are triggered instead 
by smooth density gradients — realistically computed — at 
the boundaries of the layer (top rows of Figures [3] and 

3.2. Metal-Rich Case: 4 x Solar Metallicity 

Figure [5] follows the evolution of dust that is initially 
well mixed with gas at 4x solar metallicity. It shares the 
same timeline as Figures [1] and [21 Thus the last profile 
marked / = 1.1 in Figure[S]is attained at a time 10% later 
than that marked / = 1.0 in the other figures. This last 
profile is the marginally stable state identified using our 
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Fig. 1. — Snapshots of settling dust computed with the ID 
code only. Plotted is the dust-to-gas ratio as a function of height 
at various instants of time. Relative timestamps are assigned 
by the non-dimensional parameter /; see the inset equation for 
the absolute elapsed time, which assumes p B = 1 g/cm 3 and 
F = 1 (equation [3j . Note that the shapes of these profiles and 
their relative spacing in time are independent of the absolute 
elapsed time, and thus independent of p B , F, and s. A t / = 0, 
dust begins well mixed at solar metallicity (/i = 0.015; ILoddersI 
2003). By / = 0.054, a pileup has formed at the dust layer's 
edge. The pileup is an artifact of our assumption that the 
initial dust-to-gas ratio /ii n i t is constant everywhere, which causes 
particle fluxes Pdl%ell °c lA z \ t° decrease with decreasing height \z\. 

standard procedure, for the case of supersolar metallicity. 
It achieves a midplane dust-to-gas ratio of fj,o — 20.3 — 
large enough to exceed the Toomre threshold in a disk 
that has twice the gas content of the minimum-mass solar 
nebula (F = 2 in equation [TU|) . 

The evolution of the metal-rich disk over the course of 
a total of 21 iterations — some of which are sampled in 
Figures [THIS — is similar to that of the solar metallicity 
disk, with two notable differences. When the unstably 
stratified pileups of dust collapse (iteration #4, shown in 
the top row of Figure[7]), enough dust is transferred to the 
midplane that fj, attains an appreciable maximum there. 
This bump contrasts with the nearly flat profile seen for 
the solar metallicity run (Figure [3]), and persists at least 
through iteration #16 (Figure [5]). A second difference 
is that in every KH-stable simulation following iteration 
#13, the vertical kinetic energy, although it eventually 
levels off, ends orders of magnitude higher than where it 
began (Figure[8l and top row of Figure[9]). Some currents 
and/or turbulence appear to be sustained as the state of 
marginal stability is approached. 

Related to this second point, we should acknowledge 
that our standard procedure ignores whatever velocities 
are present at the end of a given 3D simulation when ini- 
tializing the velocities of the subsequent 3D simulation. 
That is, with every iteration, velocities are set anew ac- 
cording to equation (J5]) , with vertical and radial velocities 
reset to zero. The assumption we make in our standard 
procedure is that whatever velocities are maintained in 
a KH-stable layer do not stop dust from settling at the 
local terminal velocity v rc \. A crude attempt at relaxing 
this assumption is made in ^ 14. 21 
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Fig. 2. — Snapshots of settling dust computed using the standard 
procedure of ^2]which combines the ID and 3D codes, for the case 
of bulk solar metallicity. Elapsed time is marked by /; plotted 
values coincide with those in Figure [T] The shapes of the profiles 
and their relative spacing in time do not depend on the absolute 
elapsed time; they are independent of p B , F, and s. Dust begins 
well mixed with gas at fi = 0.015 and ends in the marginally stable 
state with midplane fig = 2.45. Vertical gridpoints from the 3D 
code are plotted as dots. In comparison to the purely ID results 
of Figure [T] the pileup at the layer's edge is smoothed away, 
probably by the Ray lcigh- Taylor instability, between / = 0.22 
and / = 0.65. Except for transferring some dust at altitude 
to the midplane, the instability leaves the midplane relatively 
unaffected, which until / = 1.0 evolves much as it does in Figure[T] 

4. EXTENSIONS 

In p4.ll we test the robustness of our results against 
the size of our simulation box. In i)4.2l we modify the 
procedure of <J21 pushing to still higher dust-to-gas ra- 
tios at the midplane and revising our identification of 
marginally stable states. 

4.1. Bigger Box Runs 

Box size can artificially affect stability because a given 
box can only support modes having an integer number 
of azimuthal wavelengths inside it. Thus too small a box 
may be missing modes that would otherwise destabilize 
the layer. To assess whether our box size is too small, 
we redo the 3D simulations of our standard marginally 
stable states (iteration #19 of the solar metallicity case 
and iteration #21 of the metal-rich case), quadrupling 
simultaneously the azimuthal box size and the num- 
ber of grid points N^. By increasing both in tandem, 
we maintain the same resolution N^/L^ as that of our 
standard runs. The results are plotted as dotted lines 
in the top panels of Figures [5] and [9] For both the so- 
lar metallicity and metal-rich cases, the bigger box runs 
still yield stable layers, just as the standard box runs do. 
We conclude that our standard box sizes are probably 
adequate. 

This conclusion is a bit surprising when we compare 
our standard box size to our findings in Paper I. We 
found in Paper I that those KH modes that most visibly 
disrupted the dust layer had azimuthal wavelengths be- 
tween 2.6z max and 4.3z max . Our standard choice here for 
azimuthal box size is = 2.91z max , which at face value 
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Fig. 3. — Two successive iterations of our procedure of [[2] for 
the case of bulk solar metallicity. From left to right, the panels 
show a starting dust profile (black curve) settled by the ID code 
until its midplane fio increases by 30% (red dot-dashed curve). 
This settled curve is then passed to the 3D code and evolved (blue 
curve) until it stabilizes (rightmost panel showing how the vertical 
kinetic energy at the midplane eventually levels off). Top panels 
show iteration #6 of 19 (equivalently / = 0.33 on the timeline of 
Figure [2}- The unstably stratified pileups collapse around t ~ 11 
orbits, increasing the midplane dust content by ~20% (top middle). 
Bottom panels show iteration #7 (/ = 0.38) which begins where 
iteration #6 leaves off — except that the kinetic energy of the flow 
is reset to a low value (bottom right versus top right panels), and 
the slight asymmetry in (fi) about z = (top middle panel, blue 
curve) is dropped upon fitting a polynomial only to z > (bottom 
left, black curve). The oscillations in the /^-profiles are artifacts of 
having too few basis functions in z. They did not seem to introduce 
instability, which always occurred instead at the edges of the dust 
layer where gradients were steepest and realistically computed. 




0.01 -0.02 0.02 2 4 6 8 10 
z (H,) z (H,) time (orbits) 



Fig. 4. — Similar to Figure [3] but showing iterations #16 (top) 
and #17 (bottom) out of a total of 19, for the case of bulk solar 
metallicity. In iteration #16, dusty gas at the layer's edges mixes 
with dust-poor gas at higher altitudes (top middle), probably by 
the KHI. The subsequent evolution during iteration #17 shows no 
sign of instability after 10 orbits. 
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Fig. 5. — Similar to Figures [3] and [4] but showing the last cou- 
ple iterations (#19a and #19b) which provisionally identify the 
marginally stable state for the case of bulk solar metallicity. In- 
creasing the midplane dust content from iteration #18 by 30% 
(bottom panels) leads to a KH-unstable profile, while an increase 
of 10% preserves KH stability (top panels). Quadrupling and 
N,/, simultaneously (dotted lines) does not change our answer. The 
marginally stable state in t he top panels is refined according to a 
modified procedure in i|4.2l 
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Fig. 6. — Snapshots of settling dust computed with the full pro- 
cedure of J2]which combines the ID and 3D codes, for the case of 
4x bulk solar metallicity. Elapsed time is marked by /, measured 
on the same timeline characterizing Figures [T] and [2] The shapes 
of the profiles and their relative spacing in time do not depend on 
the absolute elapsed time; in this sense the evolution is not sen- 
sitive to p s , F, and s. Vertical gridpoints from the 3D code are 
plotted as dots. Dust begins well mixed with gas at fi = 0.06 and 
ends in the marginally stable state with midplane fio = 20.3. The 
midplane density in this last state already exceeds the threshold 
for Toomre instability in a disk with twic e the gas content of the 
minimum-mass solar nebula (equation 1101 with F = 2). 
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Fig. 7. — Two successive iterations of our procedure outlined in 
J2] for the case of 4x solar metallicity. From left to right, the 
panels show a starting dust profile (black curve) settled by the ID 
code until its midplane fio increases by 30% (red dot-dashed curve). 
This settled profile is then passed to the 3D code and evolved (blue 
curve) until it stabilizes (rightmost panel showing how the vertical 
kinetic energy at the midplane eventually levels off). Top panels 
show iteration #4 of 21 (equivalently / = 0.22 on the timeline 
of Figure [6}. When the unstably stratified pileups collapse, they 
increase the dust content of the midplane by ~35% (top middle 
panel). The resultant dust profile, settled further in iteration #5 
(bottom panels), remains free of instabilities after 10 orbits. 
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Fig. 8. — Similar to Figure[7]but showing iterations #15 and #16 
out of a total of 21 for the case of 4x solar metallicity. Shown are 
two examples of KH-stablc profiles whose midplane vertical kinetic 
energies end orders of magnitude above their starting values. Every 
3D simulation starting with iteration #13 in the metal-rich case 
shows this kind of sustained motion even though the density profiles 
may be KH stable according to our criterion. 



means that we are only resolving one wavelength at best 
of an important mode. However, this simple comparison 
may not be fair. In Paper I we studied dust layers charac- 
terized by a spatially constant Richardson number. The 
vertical density profiles there differ somewhat from those 
derived here. In particular the profiles in Paper I have 
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Fig. 9. — Similar to Figurcs[7]and[8]but showing iteration #21 in 
the top panels, in which the marginally stable state is found for the 
case of 4x solar metallicity according to our standard procedure. 
The midplane (fio) = 20.3, corresponding to a midplane density 
that exceeds that required for gravitational instability in a disk 
having twice the gas content of the minimum-mass solar nebula. 
The same dust profile inserted into a shearing box four times as 
wide in the azimuthal direction as our standard box and having 
four times as many azimuthal grid points yields qualitatively the 
same result (dotted line). Settling still further according to our 
sta ndar d procedure results in KH instability (bottom panels), but 
in i|4.2l we experiment with a modified procedure that tries to hold 
off KH instability for a while longer. 

steep density gradients everywhere, whereas here density 
gradients are steep only at the edges of the layer. When 
a layer in Paper I became KH unstable, it seemed to do 
so everywhere at once, whereas here instability always 
originates at the edges. We have verified that this is true 
even for the final iterations leading to our identification 
of the marginally stable state. Obviously these edges 
have vertical thicknesses that are smaller than that of 
the layer as a whole. Since the most unstable azimuthal 
wavelength of the KHI is expected to be of order the ver- 
tical thickness of the shearing layer fe.g.JChandra sekharl 
1981), it seems that our standard box sizes here, though 
small compared to our box sizes in Paper I, permit us 
to resolve several azimuthal wavelengths of those modes 
that are most able to disrupt the thin edges. 

4.2. Refining the Marginally Stable State Using a 
Modified Settling Procedure 

Using our standard procedure of 321 we can only pro- 
visionally identify marginally stable dust profiles. The 
identification is provisional because by settling all dust 
particles at their local terminal velocities v rc \, we wind 
up with edges so unstable that they also destabilize the 
midplane. In reality, dust particles at the edge may stop 
settling because they attain a state of marginal stability 
first, remaining lofted up by the gas motions they stir 
up locally, and leaving dust particles near the midplane 
free to keep settling. In other words, marginal stability 
may be reached sequentially, starting from the edges and 
ending at the midplane. Our standard procedure does 
not allow for this kind of gradual evolution because the 
ID code settles all dust particles at their local terminal 
speeds regardless of their location. In this sense our stan- 
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dard procedure is too blunt because it does not allow for 
slower settling at the stirred-up edges and faster settling 
at the quiescent midplane. True marginally stable states 
should have midplane dust-to-gas ratios even higher than 
the maximum ones displayed in Figures El El EJ and l9F^I 
To remedy this shortcoming, we modify our proce- 
dure by applying a weighting function < W(z) < 1 
to each dust particle's settling velocity. Starting with a 
KH-stable state near the end of our standard sequence 
of iterations, we continue to let dust particles drift to 
the midplane in the ID code, not at v TC i(z) but at the 
weighted velocity W(z)v re i(z). We use either a Fermi 
function 



W(z) 



1 



1 + exp[(z - z 50 )/z v 



(17) 



described by two parameters Z50 and z w , or a Gaussian 



W(z) = exp(-z 2 /24), 



(18) 



described by a single parameter z w . The choice of weight- 
ing function is somewhat arbitrary; it depends on the 
shape of the dust profile to be settled and is made case- 
by-case according to considerations outlined below. The 
intent of the weighting function is to slow the steepening 
of density gradients at the dust layer's edges and thereby 
prevent the edges from destabilizing the entire layer. 

We start with the KH-stable profile in iteration =#=18 
of our solar metallicity sequence (black solid curve in top 
left panel of Figure [TU)) . The dust layer is characterized 
by a "core" from z — to z « 0.5,z max over which (//,) 
is fairly constant, and an "edge" from z 0.5,z m ax to 
Z — z max over which the dust content drops to zero. Be- 
cause the instabilities that threaten to disrupt the layer 
originate in the edge and not the core, we seek a weight- 
ing function W(z) that slows the downward drift of dust 
in the edge but not in the core. At the same time W(z) 
should not be so strongly weighted toward the midplane 
that the core disconnects from the edge and opens a 
local gap in dust content. We find upon experimen- 
tation that a Gaussian does not have enough flexibil- 
ity to meet these requirements for this particular itera- 
tion. However a Fermi function with z§o — 0.005iJ g — 
corresponding approximately to the boundary of the 
core — and z w = 0.05z5o seems to work well (blue dashed 
curve in top left panel of Figure [TU]). We use this weight- 
ing function to settle the dust profile until its midplane 
fio increases by 30% to a value of 2.9 (red dot-dashed 
curve). This settled layer remains KH stable (top middle 
and right panels of Figure [TO]) — unlike the layer settled 
without the weighting function (bottom panels of Figure 

In the new profile settled with our modified procedure, 
the distinction between the core and the edge is no longer 
so sharp. Thus to settle this new profile even further, a 
simple Gaussian for the weighting function seems to suf- 
fice. Choosing z w = 0.0025_ff g ~ 0.25z max , we attempt 
to increase the midplane /to yet again by 30%, but find 
the resultant profile to be KH unstable (bottom panels 
of Figure fT0|). 

14 Another reason our dust profiles underestimate actual dust- 
to-gas ratios is because we neglect vertical self-gravi ty, which en- 
hances stability by increasing the Brunt frequency IjSekival 119981 : 
lYoudin fc Shulfeool ). 
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Fig. 10. — E xtended settling simulations using the modified pro- 
cedure of £14.21 for the case of bulk solar metallicity. We start with 
iteration #18 from our standard procedure (black curve, top left). 
A Fermi weig hting function with 250 = 0.005H g and z w = 0.05Z50 
fcquation ll7l labeled 'W at top left) allows dust near the midplane 
to settle more than dust at higher altitude. The settled profile at- 
tains a midplane (/io) = 2.9 and is KH stable (top middle and 
right panels). Further settling, this time using a Gaussian weight- 
ing function with z w = 0.0025-Hg, results in KH instability (bottom 
panels). Although the modified procedure enables us to settle be- 
yond the last stable state identified using our standard procedure, 
the gains are not large enough to reach the Toomre density in solar 
metallicity disks. 

Similar results are obtained for the metal-rich case 
(Figure QT]). Using Gaussian weighting functions we are 
able to push the midplane dust-to-gas ratio /to to a new 
record of 26.4, which is 30% greater than the value at- 
tained using our unweighted standard procedure. 

5. SUMMARY AND DISCUSSION 

To form rocky planets and gas giant cores, dust must 
amass in a circumstellar disk. In the classic scenario 
for forming planetesimals, dust settles vertically toward 
the midplane into an ever thinner and denser l ayer that 
eventually becomes grayitatio nally unstable (|Safronovl 
119691 : iGoldreich fc Wardl fl973h . Toomre's criterion for 
gravitational instability (GI) is satisfied for midplane 
dust-to-gas ratios (pd/p g )o = Mo ^ Mo, Toomre, where 
Mo, Toomre ~ 34 for a minimum-mass nebula orbiting a 
solar-mass star (equation llOl note that /io, Toomre is nearly 
independent of disk radius). For comparison, in a disk of 
well-mixed dus t and gas at solar abundance, /to ~ 0.015 
(jLoddersI 120031 ). Whether dust can accumulate until its 
density increases by more than three orders of magnitude 
depends on how turbulent the ambient gas is. Even sup- 
posing that gas in certain regions of the disk is not intrin- 
sically turbulent (e.g., because it is too weakly ionized to 
support the magnetorotational instability), the dust it- 
self may excite turbulence in gas by a Kelvin-Hclmholtz- 
type shearing instability (KHI). The KHI is triggered 
when the velocity gradient between dust-rich gas at the 
midplane and dust-poor gas at altitude becomes too 
large. Barring gravitational instability, dust should set- 
tle to a state that is marginally stable against the KHI. 
The question of whether gravitational instability is viable 
translates into the question of whether the state that is 
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Fig. 11. — Ex tended settling simulations using the modified pro- 
cedure of £14.21 for the case of 4x bulk solar metallicity. We start 
with iteration #21 from our standard procedure (black curve, top 
left). A Gaussian weighting function with z w = 0.00132H g = 
^max/12 is used to settle preferentially the innermost layers, which 
achieve a maximum (/Ml) = 26.4 and remain KH stable (top middle 
and right panels). Although further gains in fio did not material- 
ize (bottom panels, using a Gaussian of 2 W = 0.00127-Hg), fiQ is 
already high enough that gravitational instability is viable in a disk 
having ~l-2 times the gas content of the minimum-mass solar neb- 
ula. 

marginally stable to the KHI has > ^o, Toomre (this is 
a necessary but not sufficient criterion for the formation 
of planetesimals by gravitational collapse; see footnote 
E}. 

In this paper, we sought out such marginally stable 
states by numerical simulation. Starting with dust well 
mixed with gas at either bulk solar or supersolar metal- 
licity, we allowed dust to settle vertically until dynamical 
instabilities prevented the midplane density from increas- 
ing further. We tracked the approach to the marginally 
stable state using a combination of a ID settling code and 
a 3D shearing box code, working in the limit that parti- 
cles are small enough not to excite streaming instabilities. 
All the instabilities that afflicted our dust layer origi- 
nated at the layer's edges, where dust density gradients 
were steepest. We found evidence for two kinds of insta- 
bilities: the usual KHI driven by vertical shear, and the 
Rayleigh- Taylor instability (RTI) driven by the weight of 
piled- up dust. These instabilities were mostly confined 
to the dust layer's top and bottom surfaces, leaving dust 
near the midplane free to settle but occasionally speed- 
ing up the accumulation of solids by transferring dust 
from pileups downward. The midplane density stopped 
increasing when the dust layer became so thin that in- 
stabilities at the edges threatened to overturn the entire 
layer. 

Using our standard procedure of $2] we attained max- 
imum dust-to-gas ratios of /io « 2.45 and 20.3 for the 
cases of solar and 4x solar bulk metallicity, respectively 
(Figures [5] and [9]) . These values are lower limits because 
in our standard procedure dust particles at the layer's 
top and bottom faces keep settling until they excite in- 
stabilities so vigorous that dust at the midplane is stirred 
up. In reality, dust at the layer's edges may reach a state 



of marginal stability and stop settling, leaving dust near 
the midpla ne fr ee to settle further. We modified our pro- 
cedure in JO] to try to account for this effect, reaching 
Ho w 2.9 and 26.4 for the two metallicity cases (Figures 
fTOlandfTTj). These values are still lower limits because our 
simulations omit self-gravity. But the correction for self- 
gravity should be small for the solar metallicity disk, on 
the order of 10% (~ 2.9/34). For our supersolar metal- 
licity disk, the correction for self-gravity might be on the 
order of unity (~ 26.4/34)— alth ough it might also be 
much higher, as lSekival (fl998h and lYoudin fe Shul (|200l 
showed that vertical self-gravity can yield a singularity 
in hq. 

We conclude that a minimum-mass disk of bulk 
(height-integrated) solar metallicity orbiting a solar-mass 
star cannot form planetesimals by self-gravity alone: 
even neglecting turbulence intrinsic to gas, the KHI 
would force the midplane dust density to fall short of the 
Toomrc density by about an order of magnitude. Our 
results make clear what changes to the circumstcllar en- 
vironment would be needed for self-gravity to prevail. To 
attain the Toomre density in a minimum-mass gas disk, 
the bulk metallicity would need to be enhanced over solar 
by a factor of a few < 4. For disks with total mass (gas 
plus dust) greater than that of the minimum-mass solar 
nebula, the required degree of metal enrichment would 
be lower. 

Our results agre e with those of the prescriptive model 
of Weidenschilling (2006), who found that the density of 
mm-sized particles (r s ~ 0.001) at r — 3 AU in a disk for 
which p g = 1.6 x 10- 9 gcm- 3 (F w 1.3) and S d /S g w 
0.015 (solar metallicity) fell short of the Toomre density 
by about a factor of 10. When the bulk metallicity Ed/S g 
increased to 0.054, the Toomre density was exceeded by 
a factor of 3. 

5.1. How Spatially Constant is the Richardson Number? 

In Paper I as i n previous works (ISekival 119981 : 
lYoudin fc Shul l200l lYoudin fc Cmanel I2004D . all dust 
profiles were assumed to have spatially constant Richard- 
son numbers Ri. The dust profiles we have computed are 
free of this assumption, whose validity we can now test. 

We calculate Ri(z) for our marginally stable s tate s, de- 
rived under both standard ($2]) and modified f ^4.2|) pro- 
cedures. To compute the numerator (Brunt frequency) of 
Ri in equation we use the horizontally averaged dust- 
to-gas ratio (fj,(z)}, computing derivatives using centered 
differences and assuming the gas to obey a Gaussian den- 
sity profile (see footnote^. To compute the denominator 
(vertical shearing rate) of Ri, we also use (fi(z)}, insert- 
ing it into equation (J5j) and computing therefrom the 
velocity derivative. Of course we could also compute the 
denominator more directly by using the simulation out- 
put itself for — this alternative approach turns out to 
give identical results for the solar metallicity disk, but for 
the metal-rich disk the Ri{z) so generated varies much 
more erratically. As noted in fc|3.2[ the metal-rich disk 
sustains gas motions well above those we put in as seed 
noise. These motions are not strong enough to overturn 
the dust layer but they are large enough to render Ri 
highly variable, both in time and space. By not using 
the simulation data for v$(z) and relying instead on the 
better behaved (n(z)), we effectively average Ri in time 
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and space. 

Results for the solar metallicity runs are shown in Fig- 
ure [l_2j We plot Ri only where (/i)-gradients are large 
enough to be computed reliably — thus we avoid regions 
closest to the midplane. Although we find that Ri is not a 
strict constant, it varies only between 0.1 and 0.3 within 
a large fraction of the edges of the dust layer — precisely 
where instabilities, presumably shear-driven, have rear- 
ranged dust. This result compares favorably with Paper 
I, where we found that a solar metallicity disk has a crit- 
ical Richardson number of i?i cr it ~ 0.2. 

Evidence for a constant Ri within the edges of the 
dust layer is even stronger for the metal-rich disk, as 
shown in Figure [T3l Here Ri hovers near 0.5 over much 
of the edges — again consistent with Paper I. See Figure 
6 of that paper; admittedly the curve for i?i cr it(Ed/E g ) 
in Paper I needs to be extrapolated to the supersolar 
metallicity considered here, £d/E g = 0.06. 

The Ri(z) profil es in Figures [121 and [TBI differ from 
those in Figure 3 of lBai fc Stonl pOlOal) : see the dashed 
curves corresponding to their 3D simulations, all of which 
include marginally aerodynamically coupled particles. 
These differences support their arguments that their sim- 
ulations were not afflicted by the KHI. 

In summary, the assumption made in other studies 
that well-coupled particles settle into a layer for which 
Ri is spatially constant does not appear too bad. This is 
welcome news, not least because it implies that the final 
marginally stable states are robust against uncertainties 
in initial conditions ($2]). One caveat is that we have not 
tested this assumption for those regions closest to the 
midplane, as they could not relax by our method before 
being disrupted by instabilities at the edges. This is an 
area where more work can be done; see §5.31 Another 
caveat, supported independently by Paper I, is that the 
critical value of Ri to which dust relaxes is not unique, 
but increases with bulk metallicity Ed/E g . For a solar 
metallicity disk, Ri CT n « 0.2, but for one having 4x 
solar metallicity, Ri C nt ~ 0.5. This trend has not yet 
been explained. 



5.2. The Super-Linear Relation Between Maximum 
Dust-to- Gas Ratio hq 
and Bulk Metallicity Ed/E g 

The degree to which Ri is constant is related to the 
scaling between the maximum midplane dust-to-gas ra- 
tio no and bulk metallicity Ed/E g . Naively it might be 
expected that /io scales linearly with Ed/E g — a greater 
total amount of metals simply yields a proportionately 
dustier midplane — and indeed a linear relation is implied 
by our order-of-magnitude estimate in equation (JTTJ). 
But we did not find a linear trend in our simulations. 
We find instead that the relation is super-linear: a factor 
of 4 increase in Ed/E g results in a factor of 9.1 increase 
in maximum /zo (26.4 versus 2.9). 

A super-linear trend is predicted by theories as- 
suming a constant Richardson number (jSekival 119981 : 
lYoudin fc Shull20"02l ). The large gain in midplane den- 
sity afforded by a comparatively modest increase in 
bulk metallicity is partly what makes planetesimal for- 
mation by gravitational instability so attractive. In- 
creasing Ed/E g does more than just increase the total 
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Fig. 12. — Richardson numbers Ri from the marginally stable 
profile of our standard procedure (top left) and from the marginally 
stable profile of our modified procedure (bottom left), for the case 
of solar metallicity. Vertical dotted lines separate the "core" from 
the "edges" in the standard profile (top right); these dotted lines 
are extended into the bottom panels for reference. We plot Ri ev- 
erywhere except where density gradients are too small to compute 
reliably; thus we avoid the entire core region of the standard pro- 
file, and the midplane of the modified profile. Over most of the 
edges — those layers outside the dotted lines which have directly 
experienced instability, almost certainly related to the KHI — the 
Richardson number varies between ~0. 1-0.3. Thus, the traditional 
assumption that dusty sublayers relax to states of spatially con- 
stant Ri receives some empirical support from this figure. 

amount of metals in the disk — it also helps to stabilize 
it, by decreasing the vertical shear dv^/dz. In the limit 
fio ~ (£d/E g )-f/g/Az ^ 1, where Az is the dust layer 
thickness, we have from equation ([5]): 



dv<p 
dz 



T]fl K r/fJ-o riQ K r 1 



H e E d /E g 



which decreases with increasing Ed/E g . By com- 
parison the Brunt frequency [(g/ pjdp/dz] 1 ^ 2 ~ 
[(f2 K A.z / fj,) (J,/ Az] 1 / 2 ~ O k hardly changes with E d /E g . 
Thus the Richardson number increases as Ed/E g in- 
creases; the enhanced stability allows Az to decrease; 
and thus /io oc Ed/Az scales super-linearly with Ed- 

The above order-of-magnitude relations show qualita- 
tively how a super-linear trend follows from the decrease 
in vertical shear, and the consequent increase in stability, 
brought about by an increase in bulk metallicity. How- 
ever, these simple relations are not enough to quantify 
the super-linear trend because Az appears to have can- 
celled out of both dv^ j dz and the Brunt frequency. This 
difficulty is avoided in a more formal derivation of the 
relation between \io and Ed/E g , made under the assump- 
tion of constant Ri, as described in Appendix iBl 

We note further that [iq scales with the inverse of 
the radial pressure gradient parameter v max /c s (cquiv- 
alently rf) in the same super- linear way as for Ed/E g . 
The smaller is v max /c s , the greater is the maximum /iq 
attainable; the relation between these quantities is de- 
rived under the assumption of constant Ri in Appendix 
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Fig. 13. — Same as Figure IT2l except for the case of 4x bulk 
solar metallicity. Here the evidence that layers relax to states of 
spatially constant Ri is even stronger than for the case of solar 
metallicity. Moreover, the value to which Ri tends in this metal- 
rich case is higher than f or t he solar metallicity case: 0.4-0.6 here, 
versus 0.1-0.3 in Figure fl2l This trend of increasing Ri with in- 
creasing bulk metallicity S^/Sg is the same trend independently 
identified in Paper I (see Figure 6 of that paper). In the bottom 
panels showing the marginally stable profile identified using our 
modified procedure, the bumps near z ±0.005-ff g are probably 
artificial, a reflection of our imposed weighting function W(z). 

[Bl Thus we expect our numerical results for max /iq 
(2.9, 26.4) to depen d sensitively on our assumed value 
of Wmax/cs = 0.025. (|Bai fc Stond (|2010d) also reported 
that the degree of clumping caused by the streaming in- 
stability increased strongly with decreasing v max /c s .) 

5.3. Future Directions 

With each iteration of our standard procedure we al- 
lowed dust particles to settle at their full terminal ve- 
locities, regardless of gas motions evinced in previous 
iterations. We tried to account for these gas motions in 
a modified procedure by reducing settling velocities at 
altitude where dust may have already attained marginal 
stability. Settling velocities were reduced by weighting 
functions chosen by eye. This modified procedure en- 
abled us to extend the settling sequence by one iteration 
but no more. Other weighting functions might allow the 
sequence to be extended further. Introducing weight- 
ing functions earlier in the sequence (rather than at the 
end of our standard procedure, as we have done), and 
increasing the midplane density by a smaller increment 
with each iteration (less than the 30% increment that we 
have adopted), would allow for a more gradual evolution 
and possibly permit the midplane to reach still greater 
densities. 

Such a program would be straightforward to pursue 
but would be subject to arbitrariness in the forms of the 
weighting functions. A more direct approach would be to 
abandon our hybrid 1D+3D scheme and upgrade the 3D 
code to allow for a non-zero aerodynamic stopping time 
istop for dust. Then both settling and stability could be 
tracked within a single 3D code. Similar codes have been 



written (e.g., Uohansen et"al|[200l iBai fe Stond"l2010blh 
but their application has been focussed on the stream- 
ing instability, on particles having r^K^stop 0.1 and 
(model-dependent) sizes upwards of decimeters. By con- 
trast, we arc interested in the possibility that even the 
smallest particles, for which < rixistop *C 1, undergo 
gravitational instability. The problem of settling small 
particles may be coupled to the problem of settling big 
ones. Even if marginally coupled particles comprise only 
a minority of the disk's solid mass, the turbulence they 
induce by the streaming instability may prevent smaller 
particles from settling i nto the thin layers required for 
gravitational instability (|Bai fc Stondl2010al) . Quantify- 
ing what is meant by "minority" remains a forefront is- 
sue. An efficient scheme for numerically simulating this 
problem would combine the anelastic methods we have 
adopted (so that the code timestep is not limited by the 
sound-crossing time) with an implicit particle integrator 
like the kind devised by IBai fe Stond (|2010bD (so that 
the code timestep is not limited by t stop ). 

Adding self-gravity would be another improvement. As 
noted at the beginning of 33 vertical self-gravity is ex- 
pected to increase the maximum dust-to-gas ratio by of 
order 10% for the case of bulk solar metallicity. For the 
case of a few x supersolar metallicity, the magnitude 
of the correction is uncertain. It is probably at least 
of order unity, but might be much more, given the ap- 
pearance of an infin ite density cusp in t hose solutions 
of lSekival (|1998f ) and lYoudin fe Shul (|2002t) that account 
for vertical self-gravity. At the same time, self-gravity 
might actually limit maximum dust-to-gas ratios if the 
fluid becomes gravito-turbule nt without producing col- 
lapsed objects {Gammic 200l|) . 

5.4. Connection to Observations and The Need For 
Supersolar Metallicity 

Observations have unveiled several trends between stel- 
lar properties and the likelihood of planet occurrence. 
Among the most well-known is the positive correlation 
between the occurrence ra te of giant planets and the 
host s tar metallicity [Fe/Hj (|Gonz alcz 1997; S antos et al.l 
[2001 iFischer fc Valentil I2005D . I Johnson et al.l fl2oToF 
provided a comprehensive analysis, using a sample of 
1266 local stars to ask whether the trend with metal- 
licity persists across all stellar masses M*. The answer 
is contained in their Figure 2. The need for supersolar 
metallicity is clear for M dwarfs (0.2 < M*/M Q < 0.7), 
where the average metallicity of the planet-hosting stars 
is [Fe/H] = 0.4. Metal-rich stars presumably once car- 
ried metal-rich disks, and so the planet-metallicity cor- 
relatio n for M dwarfs supports our resu l ts, and those of 
others (ISekivalll998l:IYoudin fc Shull2002T; iLee et al.1 120101; 
see also Uohansen et al.1 120091 : IBai fc Stond l2010af ) that 
planetesimals form much more readily in metal-rich en- 
vironments. In particular the data for M dwarfs indicate 
that a mere factor of 10 04 = 2.5 increase in metallic- 
ity above solar substantially increases the probability of 
planet occurrence. This is consistent with our finding of 
a super-linear trend between maximum dust-to-gas ratio 
and bulk metallicity ( i)5.2l and Appendix |B|) . 

However, the planet-metallicity correlation 
weakens s ystematically with increasing stel- 
lar mass ('Joh nson et al.l I2010T ) . For A stars 
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TABLE 1 



Metallicities of Extrasolar Planets IIGuillot et al,II2006I) and Solar System Gas 

Giants (IGuillotII2005|K 



Name 


-^planet 


M z a 


■^planet 


■Zplanet / Zq 


[Fe/H]* 


■^planet / Z* 




(M ffi ) 


(M e ) 


(M z /M planet ) 








HD209458 


210 


20 


0.095 


6.35 


0.02 


6.06 


OGLE-TR-56 


394 


120 


0.304 


20.3 


0.25 


11.418 


OGLE-TR-113 


429 


70 


0.163 


10.9 


0.15 


7.7 


OGLE-TR-132 


350 


105 


0.3 


20 


0.37 


8.531 


OGLE-TR-111 


168 


50 


0.297 


19.84 


0.19 


12.81 


OGLE-TR-10 


200 


10 


0.05 


3.33 


0.28 


1.75 


TrES-1 


238 


50 


0.21 


14.0 


0.06 


12.2 


HD149026 


114 


80 


0.70 


46.78 


0.36 


20.42 


HD189733 


365 


30 


0.082 


5.479 


-0.03 


5.87 


Jupiter 


318 


10-42 


0.03-0.13 


2.0-8.8 





2.0-8.8 


Saturn 


95.2 


15-30 


0.16-0.32 


11-21 





11-21 



a The metal content for each listed extrasolar planet was derived from a model of 
a planetary interior that includes an additional energy source at the planet's center 
whose power equals 0.5% of the incident stellar luminosity. 
b The solar metallicity Z e is taken to be 0.015 (jLodd crs 2003). 



(1.4 < M*/M Q < 2.0), the correlation is arguably 
not present. This calls into question the need for 
supersolar meta llicities to form plan etesimals. The 
observations of Uohnson et al.l (|2010l ) might still be 
reconciled with gravitational instability if more massive 
stars host more massive disks, although disk mass would 
have to scale with stellar mass in a faster than linear 
way to lower the threshold Toomre density (equation 
[T0|) . The possibility also remains that the observations 
are not actually a direct or sensitive probe of the theory. 
The observations concern stellar metallicity, which 
might at best correlate with the global metallicity of the 
disk, integrated over both disk height and disk radius. 
By comparison, theory concerns the local metallicity 
Ed/Sg, integrated over height but not radius. This 
local metallicity (not to be confused with the local 
dust-to-gas ratio /i) can evolve substantially from its 
global value, as a consequence of radial particle drifts 
and photoevaporation (e.g., CY10). 

Rather than look to their parent stars for evidence for 
local disk enrichment, we can look to the planets them- 
selves. If planetesimals can only form in metal-enriched 
environments, we expect that the resu l tant p lanets will 
also be metal-enriched. iGuillot et al.l ((2006) computed 
the bulk metallicities of the first nine extrasolar plan- 
ets discovered to be transiting, all of which are hot 
Jupiters. The results are listed in Table [TJ together 
with the modeled bulk metallicities of Jupiter and Sat- 
urn. All eleven are indeed metal-enriched, by factors 
ranging from 2-47 relative to the Sun, and 2-20 rela- 



tive to their host stars. One caveat behind these re- 
sults is that models of hot Jupiter interiors are subject 
to the uncertainty over the extra source of internal heat 
responsible for their unexp ectedly large radii (see, e.g., 
iBatvgin fe Stev enson 2010), who also descr ibe a promis- 
ing so lution). To inflate planetary radii, Guill ot et al.l 
(2006) included in each hot Jupiter model an additional 
source o f power equal to 0.5% of th e received stellar irra- 
diation (jGuillot fc Showman! r2002h . The bulk metallici- 
ties inferred from the models depend on the details of this 
extra energy source. Modulo this caveat, every planet is 
enriched in metals by at least a factor of ^2 above so- 
lar, which is consistent with our finding that forming 
planetesimals by gravitational instability requires metal 
enrichments of this order. 
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APPENDIX 
BACKGROUND DISK MODEL 



For numerical estimates in this paper, we adopt the standard disk model derived in the review by lChiang fc Youdinl 
(|2010l ) . The disk has surface densities 
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in gas (g) and dust (d). The dimensionless parameters F and Z IC \ = (E/E g )/0.015, typically of order unity, describe 
how much total mass the disk has relative to the minimum-mass solar nebula, and how metal-rich the disk is compared 
with a gas of solar abundances, respectively. The mini mum-mass sol ar nebula (F = 1, Z xc \ = 1) uses a condensate 
mass fraction for solar abundances of Ed/E g = 0.015 (Loddcrs 2003 }) . Values of Z IC \ > 1 correspond to supersolar 
metallicities Ed/E g > 0.015. Integrated to r — 100 AU, equation (IA1 1) yields a total disk mass of O.O3-FM . 
At the disk midplane, the gas temperature, scale height, and density are given by 

/ r \ 

t=120 (au ) K (A3) 

F g = 0.022r(^ J ) 2/7 (A4) 

/ r \ -39/14 

Pg0 -2.7x 10- 9 F(— J gcm^ 3 . (A5) 

These are adapted from IChiang fc Goldreichl (|1997t ). adjusted for a disk obeying (IA1|) — (|A2|) . orbiting a pre- main- 
sequence star of mass M* = 1M , radius i?» = 1.7i? , and temperature T* = 4350 K. 

THE SUPER-LINEAR RELATION BETWEEN MIDPLANE DUST-TO-GAS RATIO fj, AND BULK METALLICITY S d /S g 

We derive fio as a function of Ed/E g under the assumption of a constant Ri. Some evidence supporti ng a constant 
Ri was found in our simulations f H5.ip . The density profile for constant Ri is used in a number of papers (Sckiva 1998; 
lYoudin fc Shv3[200l Paper I) and we begin by repeating the result, neglecting self-gravity as we have throughout our 
paper. The dust-to-gas ratio is given by 

1/2 

- 1 (Bl) 



where 



\/{l + ^Y + {z/z d f 



Jit Umax mo\ 

Zd = tz (B2) 

is a characteristic dust height and v max = (see equations [S] and |5J) is a constant equal to the difference in 

azimuthal velocity between a strictly Keplerian flow and dust-free gas. The dust density drops to zero at 

^0(2 + ^0) 

z = ±z max = ± — z d . (B3) 

t + Mo 

A comment on equation (IB1|) . in the limit that ^> 1: except where fi is nearly constant near z -C z max / (xq and 
where it falls to zero near z = z max , the shape of /j,(z) is that of l/z. This form follows simply from the constancy 
of Ri. Because the numerator of Ri is approximately constant with fi f ij5.2|) , the denominator must be as well: 
dv^/dz ~ (v max _/ /i) I z ~ constant, which implies [i oc l/z. From this we can deduce the super-linear trend between 
Ho and E^/Eg as follows. The integral of /i with respect to z is proportional to the total surface density of dust E^. 
Because \i oc l/z, flattening off as z decreases below z ma . x /no, this integral varies as log/io- Then fj,o ex expEd, crudely. 

More formally, we have 



Ed = 2 / pddz = 2p g0 / \x dz (B4) 
Jo Jo 

where p g o is the midplane gas density, assumed constant because z max H g . The gas density profile always well 
approximates a Gaussian (see footnote [SJ , from which it follows that E g w y/27rp s oH g . Then 

I-^jT"*- (B5) 



Inserting (|B1[) into (|B5I) we have 

= log[1 + , + M y 2(2 + Mo) i /2] _ ^f + Mo) 1 / 2 (B6) 

2 ^ Eg (1 + no) 

In the limit po ^ 1, the exponential dependence of no on Ed/E g is evident. Equation (|B6I) is plotted in Figure [T4l 
with i?i = 0.25 and v ma . x /c s — 0.025. Overlaid is the same equation but with varying Ri — Ri cl n ~ 0.25(/xo/9), the 
relation we found in Paper I (see Figure 5 of that paper) . The two data points representing the maximum fio achieved 
in this paper are also plotted. The data track the variable i?i C rit(Mo) curve much better than the constant Ri curve. 

Finally note that H s / z^ oc c s /v max enters into equation ()B6[) the same way that Ed/E g does. Thus no increases 
super-linearly with c s /w max as well. This result leads us to suspect that our numerical results for fio (2-9, 26.4) depend 
sensitively on our choice for u max /c s = 0.025. In this paper we did not run simulations with different i> max /c s and so 
did not test this suspicion. 
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Fig. 14. — Super-linear trend between the midplane dust-to-gas rati o up and height-integrated metallicity E^/Eg for dust profiles 
characterized by a spatially constant Richardson number ffi cr it- Equation |B6] is plotted twice: the dashed curve uses ffi cr it = 0.25, whereas 
the solid curve varies Ri CI it according to the relation found in Paper I: i?i cr jt ~ 0.25(/io/9) 3 ' -0 (see Figure 5 of Paper I). Both curves fix 
Dmax/cs = 0.025. The maximum values of fio achieved in this paper are plotted as po ints. These data follow the variable Ri C rit{no) curve 
more closely than the constant Ri crl t curve, corroborating the evidence we found in £15.11 that ffi C rit is spatially constant but varies with 
fio (equivalently E<j/Eg). 
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